#clear environment 
rm(list = ls())

library(readr)
library(tidyr)
require(tigris) #
require(tidycensus) #
require(sp)
require(rgdal) #
require(sf) #
require(raster)
require(stargazer) #
require(readxl)
require(haven)
require(zoo) #== for na.locf
require(snow) #
require(doSNOW) #
require(parallel)
require(foreach)
require(iterators)
require(foreign)
require(abind)

# set directories
source("./Code/set_directories.R")

# Read in source receptor matrix from Solar Siting, merge with number of systems (from Solar Siting) and Output Matrices by Pollutant
library(dplyr)
source.receptor <- readRDS(file.path(DATA.SOLAR.SITING,"CEC_to_damages_by_POLL.rds"))
source.receptor = dplyr::left_join(source.receptor, read_csv(file.path(DATA.SOLAR.SITING, "NumberSystems.csv")) %>% dplyr::mutate(zip = sprintf("%05d", zipcode)) %>% dplyr::select(zip, num_systems), by="zip") %>% replace_na(list(num_systems = 0))
source.receptor <- source.receptor %>% relocate(zip, POLL, num_systems, total_damages)

co2.source.receptor <- subset(source.receptor,POLL=="CO2")
nox.source.receptor <- subset(source.receptor,POLL=="NOx")
pm25.source.receptor <- subset(source.receptor,POLL=="PM25")
so2.source.receptor <- subset(source.receptor,POLL=="SO2")

# Multiply damages columns (one for each fips) by number of systems column
vars = names(nox.source.receptor[,4:3113])
co2.damages = co2.source.receptor %>% mutate_at(vars,~.*co2.source.receptor$num_systems)
nox.damages = nox.source.receptor %>% mutate_at(vars,~.*nox.source.receptor$num_systems)
pm25.damages = pm25.source.receptor %>% mutate_at(vars,~.*pm25.source.receptor$num_systems)
so2.damages = so2.source.receptor %>% mutate_at(vars,~.*so2.source.receptor$num_systems)

total.damages.noco2 <- data.frame(cbind(nox.damages$zip,nox.damages[,4:3113]+pm25.damages[,4:3113]+so2.damages[,4:3113]))
colnames(total.damages.noco2) <- c("zip",vars)
total.damages.noco2$zip <- as.character(total.damages.noco2$zip)

# Create data to plot maps of (received) environmental benefits in counties where they are received
total.damages.fips <- colSums(total.damages.noco2[,2:3111])
total.damages.fips <- data.frame(cbind(vars,total.damages.fips))

colnames(total.damages.fips) <- c("region","value")
total.damages.fips$region <- as.character(total.damages.fips$region)
total.damages.fips$value <- as.numeric(as.character(total.damages.fips$value))
total.damages.fips <- total.damages.fips[-1,]

# Create data to plot maps of (created) environmental benefits in zips where panels are installed (created)
total.damages.zips <- total.damages.noco2[,1:2]
colnames(total.damages.zips) <- c("region","value")

# Save damage files to use in data analysis
saveRDS(total.damages.fips, file = file.path(DATA.PROCESSED, "total_damages_fips_stdsize.rds"))
saveRDS(total.damages.noco2, file = file.path(DATA.PROCESSED, "total_damages_noCO2_stdsize.rds"))
saveRDS(total.damages.zips, file = file.path(DATA.PROCESSED, "total_damages_zips_stdsize.rds"))
saveRDS(co2.damages, file = file.path(DATA.PROCESSED, "co2_damages_stdsize.rds"))
saveRDS(nox.damages, file = file.path(DATA.PROCESSED, "nox_damages_stdsize.rds"))
saveRDS(pm25.damages, file = file.path(DATA.PROCESSED, "pm25_damages_stdsize.rds"))
saveRDS(so2.damages, file = file.path(DATA.PROCESSED, "so2_damages_stdsize.rds"))
